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(57) A computerized method and system for the ra- 
diographic analysis of bone structure and risk of future 
fracture with or with out the measurement of bone mass. 
Techniques include texture analysis for use in quantitat- 
ing the bone structure and risk of future fracture. The 
texture analysis of the bone structure incorporates di- 
rectionality information, for example in terms of the an- 
gular dependence of the RMS variation and first mo- 
ment of the power spectrum of a ROI in the bony region 
of interest. The system also includes using dual energy 
imaging in order to obtain measures of both mass and 
bone structure with one exam. Specific applications are 
given for the analysis of regions within the vertebral bod- 
ies on conventional spine radiographs. Techniques in- 
clude novel features that characterize the power spec- 
trum of the bone structure and allow extraction of direc- 
tionality features with which to characterize the spatial 
distribution and thickness of the bone trabeculae. These 
features are then merged using artificial neural networks 
in order to yield a likelihood of risk of future fracture. In 
addition, a method and system is presented in which du- 
al-energy imaging techniques are used to yield meas- 
ures of both bone mass and bone structure with one low- 
dose radiographic examination; thus, making the sys- 
tem desirable for screening (for osteoporosis and risk of 
future fracture). 
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Description 

[0001] The present invention was made in part with U.S. Government support under NIH grants/contracts CA48985 
and CA47043, Army grant/contract DAMD 1 7-93-J-3021 , and American Society grant/contract FRA-390. The U.S. 
5 Government has certain rights in the invention. 

Technical Field 

[0002] The invention relates generally to a method and system for the computerized radiographic analysis of bone 
10 structure. Specific applications are given for the analysis of the bone mass and bone pattern for the assessment of 
osteoporosis and as a predictor of risk of fracture. Novel techniques involve a directional analysis of the Fourier spec- 
trum relative to many texture measures. Additional techniques include the a one-shot dual energy exposure for the 
assessment of bone mass while simultaneously obtaining an image for the texture analysis for bone structure. 

15 Background Art 

[0003] Osteoporosis is a widespread medical condition that affects about 15-20 million people in the United States 
and accounts for about 1 .3 million new fractures per year in people greater than 45 years of age. Osteoporosis man if ests 
as a loss in bone mass, a tendency to fracture and as a structural alteration of bone. Quantitative measures of bone 

20 mass serve as important diagnostic indicators for determining the risk of fracture and in following the progress of 
patients on therapy for osteoporosis. The most widely used methods of assessing bone mass are by bone mineral 
densitometry (BMD) that involves dual photon absorptiometry with either an x-ray or nuclear source, and quantitative 
computed tomography. These methods are very accurate in determining bone mass, which has been shown to be a 
very good predictor of fracture risk in osteoporosis. There is, however, considerable overlap of the measurements of 

25 BMD for patients with osteoporosis who have, or go on to have atraumatic fractures compared to age-matched controls 
who do not have, or do not go on to have, atraumatic fractures. In addition to bone mass, bone structure is probably 
also important in determining the mechanical strength of bone and thus fracture risk. A few preliminary studies have 
been performed in relating certain textural measures to bone structure. 

30 Disclosure of the Invention 

[0004] Accordingly, an object of this invention is to provide a computerized method and system for the radiographic 
analysis of bone structure and risk of future fracture. 

[0005] Another object of this invention is to provide a method and system for texture analysis for use in quantitating 
35 the bone structure and risk of future fracture. 

[0006] Another object of this invention is to provide a method and system for incorporating directionality information 
in the analysis of the bone structure (texture). 

[0007] Another object of this invention is to provide a method and system for using dual energy imaging in order to 
obtain measures of both bone mass and bone structure with one exam. 

40 [0008] These and other objects are achieved according to the invention by providing a new and improved method 
and system for the analysis of bone structure and future risk of fracture. Specific applications are given for the analysis 
of regions within the vertebral bodies on conventional spine radiographs. Techniques include novel features that char- 
acterize the power spectrum of the bone structure and allow extraction of directionality features with which to charac- 
terize the spatial distribution and thickness of the bone trabeculae. These features are then merged using artificial 

45 neural networks in order to yield a likelihood of risk of future fracture. In addition, a method and system is presented 
in which dual-energy imaging techniques are used to yield measures of both bone mass and bone structure with one 
low-dose radiographic examination; thus, making the system desirable for screening (for osteoporosis and risk of future 
fracture). 

50 Brief Description of the Drawings 

[0009] A more complete appreciation of the invention and many of the attendant advantages thereof will be readily 
obtained as the same becomes better understood by reference to the following detailed description when considered 
in connection with the accompanying drawings, wherein: 

55 

FIGURE 1 is a schematic diagram illustrating the method for analysis of bone structure according to the invention; 
FIGURE 2 is a schematic diagram illustrating the placement of ROIs on the vertebral bodies in digital lumbar spine 
images. 
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FIGURE 3 is schematic illustrating corrections for the possible nonlinear nature of the detector system's charac- 
teristic response (H & D curve for film) and for background within the ROI image data. 

FIGURE 4 is a schematic diagram illustrating the power spectrum (with sectors indicated) obtained from the Fourier 
transform of the corrected ROI image data. 
5 FIGURE 5 is a schematic diagram listing the various measures including directionality measures obtained from 

the power spectrum of the image data. 

FIGURE 6 is a schematic illustrating some of the texture measures for nonosteoporotic ("healthy") bone and for 
diseased bone. 

FIGURE 7 is a graph showing measures of bone mass for 43 patients: some with a fracture elsewhere in the spine 

10 and some without fracture. 

FIGURE 8 is a graph showing the relationship between BMD measures (bone mass) and RMS variation (bone 
structure) for patients: some with a fracture elsewhere in the spine and some without fracture. 
FIGU RE 9 is a graph showing the relationship between BMD measures (bone mass) and first moment of the power 
spectrum (bone structure) for patients: some with a fracture elsewhere in the spine and some without fracture. 

15 FIGU RE 1 0 is a graph illustrating the relationship between RMS variation and first moment of the power spectrum 

for ROIs selected from the L3 level for patients with and without fracture elsewhere in the spine. 
FIGURE 11 is a graph illustrating the relationship between the standard deviation of the angular dependence of 
the RMS variation and the minimum value of the angular dependence of the first moment of the power spectrum 
for ROIs selected from the L3 level for patients with and without fracture elsewhere in the spine. 

20 FIGURE 12 is a graph showing ROC curves calculated for the measures of bone mass (BMD), RMS variation and 

first moment of the power spectrum. 

FIGURE 13 is a graph showing ROC curves calculated for the measures of bone mass (BMD), the standard de- 
viation of the angular dependence of the RMS variation and the minimum value of the angular dependence of the 
first moment of the power spectrum. 
25 FIGURE 14 is a graph showing the average values for the texture measures for cases with fracture elsewhere in 

the spine and for cases without fracture. 

FIGURE 15 is a graph indicating the performance of the individual texture measures in the task of distinguishing 
those cases with fractures elsewhere in the spine from those without fracture. 

FIGURE 16 is schematic diagram of the artificial neural network used in merging the various bone structure features 
30 into a likelihood of risk of future fracture. 

FIGURE 17 is a graph showing ROC curves calculated for three neural network combinations. Two of the combi- 
nations include measures of both bone mass and bone structure; one of the combinations includes only measures 
of bone structure. 

FIGURES 1 8A und B show two graphs indicating the histogram (distribution) of the output values from the artificial 
35 neural network for two of the neural network combinations. 

FIGURE 19 is a schematic block diagram illustrating a system for implementing the method for the computerized, 
radiographic analysis of bone structure and risk of future fracture. 

FIGURES 20 A to D contain tables showing the effect of pixel size on four of the texture measures in terms of Az 
in predicting fracture elsewhere in the spine. 
40 FIGURE 21 is a schematic illustrating a method for the measurement of both bone density (bone mass) and bone 

structure from a single-projection, dual-energy radiographic image of the some bony body part such as the spine, 
hip or extremities according to the invention. 

FIGURE 22 is a schematic diagram illustrating two possible methods of obtaining the dual-energy radiographic 
images. 

FIGURES 23A and B are a schematic diagrams illustrating two possible methods for energy subtraction as it relates 
to the measures of bone mass. 

FIGURE 24 is a schematic diagram illustrating one possible calibration method for measuring bone mass from 
dual-energy images, including compensation for the heel effect and calibration for body thickness. 
FIGURE 25 is a schematic diagram illustrating a calibration phantom and its positioning during a patient exam. 
so FIGURES 26A to D are graphs showing the relationship between gray values obtained from low and high energy 

images of a lucite phantom; illustrating the calibration method. 

FIGURE 27 is a graph showing the relationship between the measured values for bone mass (from a BMD phantom) 
and the accepted values for the particular phantom, and thus indicating the potential for using the technique for 
measuring bone mass (along with bone structure). 
55 FIGURE 28 is a schematic block diagram illustrating a system for implementing the method for the radiographic 

analysis of both bone mass and bone structure, and thus, risk of future fracture. 

Figure 29 illustrates the logarithmic relationship between "surface area" and effective pixel length for a ROI. Two 
distinct linear portions are demonstrated and each is fitted with a straight line of different slope, which is used to 
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calculate fractal dimension. 
Best Mode for Carrying Out the Invention 

[0010] Referring now to the drawings, and more particularly to Figure 1 thereof, a schematic diagram of the analysis 
of bone structure is shown. In this example, the aim is to extract the characteristics of the bone trabeculae using texture 
analysis of image data from digital images of bony parts of the body such as the spine. The overall scheme includes 
an initial acquisition of a radiograph of the spine (step 10) and digitization (step 20) (or a direct digital acquisition of 
the radiographic image of the spine). A region of interest (ROI) is then placed over a vertebral body on the image and 
the corresponding image data are stored in memory (step 30). Background trend correction (step 40) is performed to 
yield the underlying fluctuations, i.e., the trabecular pattern. The image data in the ROI are then input to a Fast Fourier 
Transform (step 50) and the power spectrum is calculated (step 60). Various textures measures are calculated from 
the power spectrum data (step 70) and these are merged using an artificial neural network (step 80) to yield a likelihood 
of risk of future fracture (step 80). Other texture analyses can be used such as fractal analysis. 
[0011] Figure 2 illustrates the placement of ROIs on the vertebral bodies in the digital lumbar spine images. Shown 
here are ROIs, 64 pixels by 64 pixels in size, placed at the L2, L3, and L4 levels on the spine. Placement is performed 
such that the ROIs avoid overlapping edges, bowel gas, and soft tissue folds. In general, ROIs placed at the L3 level 
had the least interference from edges and bowel gas, and thus precise placement of the ROIs within the vertebral body 
is not necessary at the L3 level. 

[0012] Figure 3 illustrates the corrections for the possible nonlinear nature of the detector's characteristic response 
(the H & D curve for radiographic films as detector) and for the background trend within the ROI image data. Background 
trend connection is necessary since the variation in optical density within the ROI in spine images includes that due 
to the gross anatomy of the human body and to the presence of bowel gas (background treads) and that due to the 
fine underlying texture which is related to the trabecular pattern of the bone. The nonuniform background trend can 
be determined using a 2-dimensional surface fitting technique (such as one with a second degree polynomial function). 
The fitted trend is subtracted from each ROI in order to yield the underlying fluctuations, i.e., the trabecular pattern. 
[0013] Figure 4 illustrates the power spectrum of ROI image data. The axes are in terms of spatial frequencies. It 
should be noted that strictly speaking, however, the power spectrum needs to be determined from an ensemble average 
of the square of the Fourier transform over an infinitely large area. The sectors indicate the method used in dividing 
the power spectrum into pie-shaped sections. Texture measures are calculated for the entire power spectrum as well 
as for the individual sectors, thus, yielding directionality measures. The power spectra of the trabecular bone pattern 
may contain low-frequency components due to some residual uncorrected background trend and very high-frequency 
components due to radiographic mottle in the original bone radiographic image. Thus, the power spectra may be filtered 
by the human visual system response function with acts as a band-pass filter. 

[0014] Figure 5 is a schematic diagram listing the various measures obtained from the power spectrum data. The 
texture analysis process initially involves two measures: the root-mean-square (RMS) variation (I RMS or shorthand, 
R) and the first moment of the filtered power spectrum (IFMP or shorthand, M), which represent the magnitude and 
the coarseness of trabecular pattern, respectively. These measures are given by 



where V(u,v) and F(u,v) correspond to the visual system response and the Fourier transform of the trabecular pattern, 
respectively. Higher moments of the power spectra can also be calculated. Higher moments are not conceptualized 
visually as easily as the RMS variation and first moment values, however. 

[0015] Due to the strong directional appearance of trabecular patterns, the RMS variation and first moment of the 
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power spectra will be calculated also as a function of angle in the Fourier domain as given below by the inequalities 
and tables. 

V©i * 9 < © 2 ) =/£ S KJ 2 for 0 L s tan-M-) < ® 2 

V m n /H 

[0016] Angular dependence of RMS variation: 

[0017] Angular dependence of First Moment of the Power Spectrum: 
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[0018] The angular dependence of the two measures (IRMS and IFMP) is examined by dividing the power spectrum 
into several sectors and performing the summations within each sector. From studies, we have found that those with 
fracture elsewhere in the spine exhibit a higher minimum value of the angular dependence of IFMP. By taking the 

20 minimum value we force the directionality measure (i.e., perpendicular to the trabeculae) for normal patients since the 
bone trabeculae without osteoporosis is assumed to not show a "washed-out" appearance and thus the directionality 
is strong as schematically shown in Figures 6 A and B. Since the trabeculae are not "washed-out" for normal patients, 
their spatial distribution would contain lower frequency structures in a direction perpendicular to the trabeculae. Oste- 
oporotic patients would tend to exhibit a more isotropic distribution due to the washed-out appearance of the trabeculae. 

25 Edge gradient analysis on the ROI data can also be performed to extract the direction perpendicular to the major 
trabeculae. The angle exhibiting the largest cumulative edge gradient is expected to indicate the direction perpendicular 
to the major trabeculae within the ROI. In addition, due to the possibility that quantum mottle and x-ray scatter may 
"hide" the underlying texture pattern of the bone trabeculae, the power spectra of uniform tissue regions within the 
medical image are also determined and used to normalize the power spectra obtained from the ROIs in the bony 

30 regions, prior to calculation of the texture measures. These analyses are expected to be useful in analyzing both the 
primary (approximately horizontal) trabeculae and the secondary (approximately horizontal) trabeculae and the sec- 
ondary (approximately vertical) trabeculae. 

[0019] Studies were done using 43 patient cases in which some had a fracture elsewhere in the spine and some did 
not. This method for evaluation was used since the texture measures here are being examined at one point in time 

35 and it has been shown that the presence of pre-existing vertebral body fractures is a powerful predictor of future risk 
for vertebral body fracture. Figure 7 is a graph showing the distribution of BMD measures (bone mass) for patients 
with fracture elsewhere in the spine and for those without fracture. Notice that the BMD values are low in nearly all the 
fracture cases as expected; most of the nonfracture cases also, however, have low BMD values, thus, demonstrating 
the need for a measure with higher specificity. 

40 [0020] Figures 8 and 9 demonstrate the relationships between BMD measures (bone mass) and the RMS variation 
for the same patients and between BMD and the first moment of the power spectrum. It is apparent that there is not a 
strong correlation between bone mass and bone structure using at least these measures for bone mass and bone 
structure. Note that in the following example, normalization of the power spectra was not included. 
[0021 ] Figure 1 0 is a graph illustrating the relationship between the RMS variation and the first moment of the power 

45 spectrum for ROIs selected from the L3 level for patients with and without fractures elsewhere in the spine. It is apparent 
that patients with fractures elsewhere in the spine tend to have a high first moment measure and a low RMS variation. 
Figure 11 is a graph illustrating the relationship between the standard deviation of the angular dependence of the RMS 
variation and the minimum value of the angular dependence of the first moment of the power spectrum for ROIs selected 
from the L3 level for patients with and without fractures elsewhere in the spine. 

50 [0022] Figure 12 is a graph showing ROC curves calculated for the measures of bone mass (BMD), RMS variation 
and first moment of the power spectrum. Here the ROC analysis was performed with respect to the task of determining 
whether or not the patient had a fracture elsewhere in the spine. The Az values (area under the ROC curve) for RMS 
variation and the first moment are superior when compared to the Az value for the measure of bone mass (BMD). 
Figure 13 is a graph showing ROC curves calculated for the measures of bone mass (BMD), the standard deviation 

55 of the angular dependence of the RMS variation and the minimum value of the angular dependence of the first moment 
of the power spectrum. 

[0023] Figure 1 4 is a graph showing the average values for the texture measures for cases with fracture elsewhere 
in the spine and for cases without fracture. Note that the values have been normalized between 0 and 1 . Figure 15 is 
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a graph indicating the performance of the individual texture measures in the task of distinguishing those cases with 
fracture elsewhere in the spine from those without fracture. Note the higher the Az value the better the performance. 
[0024] Once the texture measures are calculated they can be merged using an artificial neural network (ANN) in 
order to yield a likelihood of future risk of fracture. Figure 16 is a schematic diagram of the artificial neural network 
5 used in merging the various bone structure features into a likelihood of risk of future fracture. 

[0025] Figure 1 7 is a graph showing ROC curves calculated for three neural network combinations. Two of the com- 
binations include measures of both bone mass and bone structure; one of the combinations includes only measures 
of bone structure. 

[0026] Figures 18A und B shows two graphs indicating the histogram (distribution) of the output values from the 
10 artificial neural network for two of the neural network combinations. The output from the neural network can bethresh- 
olded so that only cases with a certain value from the neural network output are noted as having a higher risk for future 
fracture. 

[0027] Figure 19 is a more detailed schematic block diagram illustrating a system for implementing the method of 
the invention for analysis of the bone trabecular structure. Referring to Figure 1 9, radiographic images of an object are 

15 obtained from an image acquisition device (1000) and input to the system. Each bone image is digitized and put into 
memory (1001). If the radiographic image is obtained with a direct digital device then there is no need for digitization. 
The image data is first passed through the ROI selection circuit (1 002), the nonlinear detection system correction circuit 
(1003) and the background trend correction circuit (1004). The data is passed to the FFT circuit (1005) and the power 
spectrum circuit (1006). Power spectrum data are passed to the texture measure circuit (1007) and to the optional 

20 ANN circuit (1 008) in order to determine the likelihood for risk of future fracture, during which time the data are retained 
in image memory (1009). In the superimposing circuit (1010) the results are either superimposed onto images, stored 
in file format, or given in text format The results are then displayed on the display system (1020) after passing through 
a digital-to-analog converter (1030). 

[0028] The particular image acquisition system used may affect the texture measures, so the ability of the compu- 
25 terized scheme to assess osteoporosis and risk of fracture as a function of pixel size of the digitization system was 
investigated. Use of a smaller pixel size allows higherspatial frequency components to be examined. The results shown 
earlier were obtained from digitization of film at a pixel of 0.175 mm with 10-bit quantization. If the texture measures 
can still be reliably obtained at larger pixel size then direct digital systems for imaging the bone will be more easily 
produced. Figures 20A to D contains table showing the effect of pixel size on four of the texture measures in terms of 
30 Az in predicting fracture elsewhere in the spine. Results are given for variations in the aperture size (blur) and sampling 
distance for the same 43 cases used in the earlier samples. 

[0029] Figure 21 is a schematic illustrating a method for the measurement of both bone density (bone mass) and 
bone structure from a single-projection, dual-energy radiographic image of some bony body part such as the spine, 
hip, or extremities according to the invention. Such a system produces a high-energy image and a low-energy image 

35 from either a "one shot" exposure technique that employs two detectors sandwiched together or a "two- exposure" 
technique that utilizes two exposures to the patient. Figures 22A and B are schematic diagrams illustrating two possible 
methods of obtaining the dual energy radiographic images, Such systems utilize "energy subtraction" techniques to 
yield "bone-cancelled" images and "soft-tissue -cancel led" images. Such "dual-energy" systems have been developed 
for moderately-high-spatial-resolution imaging of the soft tissue in chest and for very-low-spatial-resolution acquisition 

40 of bone mass (BMD) in, for example, the spine. However, such "dual-energy" systems have not been developed for 
moderately-high-spatial-resolution imaging of bone due to the large amount of quantum mottle that results in the bone 
image (i.e., the soft-tissue-cancelled image). Moderately high spatial resolution is desirable forthe analysis of the bone 
structure, though in the past, only bone mass was of interest, and thus, the low-resolution system was acceptable. 
However, now with the new method presented earlier in this invention application that yields measures of bone structure, 

45 it is desirable to have a system that can measure both bone mass and bone structure (as opposed to subjecting the 
patient to two examinations: one for bone mass (BMD) and one for bone structure (spine radiograph).) The following 
presents such a system using computed radiography as the detector in the example. However, the method is not limited 
to computed radiography as the detector. 

[0030] Computed radiography (CR) is a digital radiographic imaging system that uses plates consisting of BaFBr 
so phosphor, a stimulable phosphor, to image the radiographic x-ray image. The pixel value in a CR image can be con- 
verted directly into x-ray exposure. The method uses dual-energy computed radiography (CR) imaging to obtain bone 
mass in a manner quite similar to that of DXA (BMD). Differences include a dependence on scatter due to the fact that 
conventional radiographic spine images are obtained with a broad area beam whereas the DXA scans are obtained 
with a low-resolution, pencil-beam geometry. However, the CR images are of high spatial resolution thus allowing for 
55 the low-energy image to be used for the texture measures of bone structure. Note that texture analysis is not possible 
on the tissue-cancelled images due to the presence of large quantum mottle (and the inability of increasing the exposure 
a lot due to patient dose considerations). The measure for bone mass is performed in a way that the region on the 
bone image that encloses the spongiosa will be integrated. All this accomplished with just one examination. 
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[0031 ] In this example, dual-energy bone images of the spine, hip and extremities are obtained using the CR system 
and the "one-shot" exposure technique. The method uses conventional x-ray equipment to produce "bone" and "soft 
tissue" images in exact temporal and spatial registration with a single x-ray exposure. Use of the one-shot technique 
eliminates motion artifacts between the high- and low-energy images and also avoids rapid switching of the x-ray tube 
5 voltage. 

[0032] With the one-shot technique, it is advantageous for the input x-ray spectrum to be double peaked. Thus, K- 
edge filtration is used to produce a double-peaked x-ray spectrum. The CR plates consist of BaFBr phosphor, and thus 
the broad beam x-ray spectrum emitted from the x-ray tube is prefiltered so that the absorbed spectrum for the front 
CR plate will peak in the region of high barium attenuation coefficient. A prefilter of 300 mg/crrV gadolinium can be 
10 used. In order to compensate for the attenuation of the x-ray tube output by the K-edge prefilter, the tube loading (mAs) 
of the x-ray tube is increased, 

[0033] Another filter is sandwiched between two imaging plates (made of the stimulable phosphor) having wide 
exposurelatitude characteristics. The front CR plate of the sandwich preferentially absorbs low-energy x-ray photons 
and transmits high-energy photons. The high-energy photons are absorbed partially in the back plate yielding two 
15 simultaneously acquired images with different effective energies. The filter serves to increase the effective energy of 
the x-ray spectrum incident on the second imaging plate. Readily available materials for this filter are copper foil or the 
CR plates themselves (which contain barium). In the results presented later in this application, the filter consisted of 
two CR plates (200 mg/cm2 BaFBr) for this intermediate filtration. 

[0034] Figures 23A and B are schematic diagrams illustrating two possible methods for energy subtraction as it 

20 relates to the measures of bone mass. In method A, the low-energy image and the high-energy image are first regis- 
tered, passed through energy subtraction and then ROIs in the bone image that are within the vertebral body are 
integrated to yield a measure related to bone mass. In method B, the low-energy image and the high-energy image 
are each separately subjected to ROI placement and integration, and then a weighted sum of the two integrated values 
(with corrections for patient body size, scatter radiation present, etc.) is calculated to yield a measure related to bone 

25 mass. An advantage of method B is that image registration in the conventional (high resolution) sense is not necessary. 
However, a way (such as manual placement) to ensure location of corresponding lumbar ROIs on the low and high 
energy images is necessary. In this application are presented results using bone phantoms. The pair of digital images, 
obtained from the two CR plates in the sandwich cassette, are read digitally by the CR system. 
[0035] Figure 24 is a schematic diagram illustrating one possible calibration method for measuring bone mass from 

30 dualenergy images, including compensation for the heel effect and calibration for body thickness. Figure 25 is a sche- 
matic diagram illustrating a calibration phantom and its positioning during a patient exam. In the example, two phantoms 
were used. One of these phantoms is used to calibrate and consists of three cylinders of synthetic bone material. The 
other phantom was also made of synthetic bone, but was shaped to look like four lumbar vertebrae and encased in 
lucite. Lucite was added to the top of the phantoms to simulated additional soft tissue, i.e., patients of varying thickness. 

35 The phantoms were imaged with the one-shot, dual-energy technique and quantities such as thickness of lucite and 
energy of the x-ray beam were varied in different trials. 

[0036] Figures 26A to D show graphs of the relationship between gray values obtained from the low and high energy 
images of the phantom. The linear fit of these data is used to determine the weights for the weighted sum of the 
integrated value of the ROI on the low-energy image with that of the ROI on the high-energy image. Values from these 
40 graphs are then used in a look up table for different imaging techniques (i.e., kVp) and for patients of different thick- 
nesses. The value obtained from the weighted sum is related to the bone mass. In a particular region of the image, 
both bone and soft tissue contribute to the attenuation of the x-ray beam. The amount proportional to the thickness of 
the bone can be determined by a weighted sum (or could be thought of as a weighted subtraction), pixel by pixel, of 
the high and low energy image data, namely, 

45 

B(x, y) = L(x,y)-W*H (x, y), 

where x,y is the location of the pixel, and L and H are the values in the low and high energy images, respectively. W 
so js the weight determined from the slope of the linear fit, as demonstrated in Figure 26. For bone mass, the integration 
of a region on a bone-only image is of interest. Thus, the method either does the weighted summation on the image 
data as shown in the above equation and then integrates on the noisy bone image (in which it may be difficult to define 
edges of the bone, such as the edge of the vertebral body) or integrates on the low-energy and high-energy images 
separately prior to the weighted summation, namely 

55 

B = sum over ROI of [L(x,y)] - W • sum over ROI of [H (x, y)] t 
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where B here is proportional to the bone mass within the ROI. The edges of the vertebral bodies are easier delineated 
on both the low-energy and high-energy images, thus making locating the ROI easier. 

[0037] Figure 27 is a graph showing the relationship between the measured values for bone mass (from a BMD 
phantom) and the accepted values for the particular phantom, and thus indicating the potential for using the technique 

5 for measuring bone mass. Since the measured values are not scaled, one can only look at the general trend of the 
data. The measured values from the weighted sum of the integrated ROIs from the low- and high-energy images follow 
the same order as does the BMD measures obtained from a Lunar DPX system. This is especially so in the images 
with only 10 cm of lucite on top of the phantom. Those with 20 cm were less stable, as expected due to increased 
scattering. This can be improved with the use of better antiscatter grids in the radiographic protocol or with a direct 

10 digital acquisition device (which have been shown to be 99% in scatter rejection). 

[0038] The computerized texture analysis for bone structure (that was presented earlier in this application) is then 
performed on the low-energy images in order to measure the quality and architecture of the bone trabeculae. The 
texture measures are not determined on the bone images (i.e., tissuecancelled images) since the large amount of 
radiographic mottle can "hide" the underlying texture of the bone structure. Thus, this dual-energy technique allows 

is quantization of both bone mass and bone structure, as demonstrated earlier in Figure 21 . 

[0039] Figure 28 is a more detailed schematic block diagram illustrating a system for implementing the method of 
the invention for analysis of both the bone mass and the bone trabecular structure. Referring to Figure 28, two radio- 
graphic images (low-energy and high-energy) of an object are obtained from an image acquisition device 2000 and 
input to the system 2100. Each bone image is digitized and put into memory (2001). If the radiographic images are 

20 obtained with a direct digital device then there is no need for digitization. The image data are first passed through the 
ROI selection circuit (2002), the nonlinear detection system correction circuit (2003), and is then passed to the entry 
circuit for bone mass (2004) and to the entry circuit for bone structure (2005). For bone mass, the data are passed to 
the integration circuit (2006) and calibration circuit (2007). From there, the data are passed to the weighted sum circuit 
(2008) and saved in memory (2009). For bone structure, the data from the lowenergy image are passed from 2005 to 

25 the background trend correction circuit (201 0). The data is passed to the FFT circuit (2020) and the power spectrum 
circuit (2030). Power spectrum data are passed to the texture measure circuit (2040). Measures of both bone mass 
and bone structure are then passed to the ANN circuit (2050) in order to determine the likelihood for risk of future 
fracture, during which time the data are retained in image memory (2060). In the superimposing circuit (2070) the 
results are either superimposed onto images, stored in the file format, or given in text format. The results are then 

30 displayed on the display system (2080) after passing through a digital-to-analog converter (2090). 

[0040] Fractal analysis has also been preliminarily investigated as a means of analyzing the bone texture. Using the 
bone ROIs shown in Fig 2, the surface area of each of the ROIs was computed at 7 different levels of resolution using 
pixel gray level analogously to height in the fractal calculations. Different levels of resolution (pixel length) were obtained 
by successively averaging larger and larger groups of adjacent pixels. Figure 29 shows the log plot of surface area vs. 

35 pixel length for one ROI. The presence of the two distinct linear portions suggests a multifractal structure. Slopes of 
the overall graph of each ROI, as well as the slopes of each of the two portions of each graph were then used to obtain 
an estimate of the overall fractal dimension as well as an estimate of the fractal dimension at stronger and weaker 
levels of resolution corresponding to the 2 distinct portions of the graphs. Using ROC analysis with the fractal dimension 
of each ROI as the decision variable, the Az obtained using overall fractal dimension was 0.65, using the fractal di- 

^o mension at the finer resolution portion was 0.76 and using the fractal dimension at the coarser resolution was 0.87 as 
compared to an Az of 0.6 obtained using BMD. 

[0041] Multifractional analysis can also be used to characterize the bone texture within the ROIs. The fractal dimen- 
sion of these ROIs will be estimated using a surface area technique, modified from one described for the computerized 
analysis of mammograms. The gray level of each pixel will be regarded as a "height" with pixel size as "length" and 
45 "width" to calculate a "surface area" for each ROI. Adjacent pixels will be then combined to yield an effectively larger 
pixel size with a new gray level averaged from these combined pixels. A new "surface area" will be calculated for each 
ROI, and the process will be successively repeated, combining adjacent pixels from earlier steps, and calculating the 
resultant surface area for each new effective pixel size. The fractal dimension (D) for each ROI is calculated, using: 



where H is the slope of a least-squares line fitted to a plot of log surface area versus log pixel size for each ROI. The 
number 2 is the topological dimension of the gray level surface. The plot (as we have found) may exhibit a multifractal 
55 nature by indicating two linear regions - a textural (fine) fractal dimension and a structural (coarser) fractal dimension. 
Both the fine and the coarse fractal dimensions can be used as texture measures. 

[0042] In two preliminary studies using separately the ROIs of the spine described above and ROIs from normal and 
osteoporotic hands, artificial neural networks (ANN) were employed. The input to the neural network was the normalized 
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power spectrum data from the background-corrected ROI. Using ROIs of size 32 by 32 pixels, the successful ANN 
contained 512 (32 * 16) input units, 40 hidden units and one output unit. The value of the output unit served as the 
decision variable. The ANN was trained using an output of 1 for abnormal (osteoporotic) and 0 for normal. Using the 
43 cases, the ANN successfully classified all abnormal ROIs as osteoporotic. 

5 [0043] Artificial neural networks (ANN) can also be applied to the differentiation of texture patterns of bone trabeculae. 
ANN is a non-algorithmic approach to information processing. Unlike many artificial intelligence techniques, which 
require extensive knowledge of the many parameters involved, ANNs learn directly from examples that are provided 
repeatedly. Once trained, a neural network can distinguish among input patterns on the basis of its learning experience. 
The analysis of the texture patterns will be performed using the image data in the spatial frequency domain in order 

10 to eliminate the shift-variant nature of the image data. The ROIs will be evaluated by calculating the power spectra by 
Fourier transformation of the background-corrected ROIs and scaled. The scaled power spectra will then be used as 
input to the neural network. Thus, for ROIs of size 32 by 32 pixels, the resulting number of input units is 1 6 by 32; due 
to the symmetry in the Fourier transformation and subsequent calculation of the power spectrum. A three-layer, feed- 
forward neural network can be used with one output unit. A back-propagation algorithm with generalized delta rule will 

15 be employed in the training process. The input, which corresponds to the corresponding power spectrum, is provided 
to the input layer of the neural network. The desired output (truth) is provided to the output layer during the training 
process. The hidden layer, which is the key element in mapping of the input patterns to the output values, is located 
between the input and output layers. A nonlinear logistic function will be used as the activation function for each process- 
ing unit in the neural network, in which 

20 

O 1 

PJ 1 + expCLfWjP^Qj) 

25 where O pj is the jth element of the actual output pattern produced by the presentation of input pattern p, Wj, is the weight 
from the ith to the jth units, and 0j is the threshold of the jth units. In the training process, the internal parameters of 
the connections between layers (including threshold values of each unit) are adjusted iteratively so that the difference 
between the output values and the desired results is minimized. This can be accomplished by the following rule: 

30 

AWj,(n+1) = Tl(SpPpi) + aAWji(n), 

where n indexes the number of iterations, -n, is the learning rate, 5 p | is the error signal, which is related to the difference 
between the output of the neural network and the desired output, and a is a momentum term that determines the effect 
35 of past weight changes on the current direction of movement in weight space. The desired "truth" for use in training 
the ANN will initially be either a 1 or a 0, where 1 corresponds to the patient having a fracture elsewhere in the spine 
and 0 corresponding to the patient not having such a fracture. 

[0044] Obviously, numerous modifications and variations of the present invention are possible in light of the above 
technique. It is therefore to be understood that within the scope of the appended claims the invention may be practiced 
^o otherwise than as specifically described herein. Although the current application is focused on radiographic medicinal 
images, the concept can be expanded to segmentation in other images of the human body. 



Claims 

45 

1 . A method for analyzing bone, comprising: 

obtaining an image containing said bone from a memory; 
selecting a region of interest of said bone; 
so determining at least one trabecular texture measure of said region of interest of said bone; and 

analyzing said bone using said at least one trabecular texture measure, wherein determining said at least one 
texture measure comprises: 

determining a power spectrum of said region of interest; and determining a root mean square variation of 
55 said power spectrum and determining an angular dependence of said variation ; and/or 

determining a first moment of said power spectrum and determining an angular dependence of said first 
moment. 
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A method as recited in Claim 1 , comprising: 

background-trend correcting said region of interest to produce a background-trend corrected region of interest; 



determining said at least one texture measure of said background-trend corrected region of interest. 
A method as recited in Claim 1 , wherein selecting said region of interest comprises: 

selecting said region of interest on at least one vertebral body in a spine. 
A method as recited in Claim 1 , wherein determining said at least one texture measure comprises: 

determining said at least one texture measure using said power spectrum. 

A method as recited in Claim 4, comprising: 

dividing said power spectrum into a number of sectors; and 
determining said at least one texture measure for each of said sectors. 

A method as recited in Claim 4, wherein determining said at least one texture measure comprises: 

determining a root mean square variation and a first moment of said power spectrum. 

A method as recited in Claim 6, comprising: 

filtering said power spectra; 

determining said root mean square variation as R: 



and 




determining said first moment as M: 




_ - -oo-oo 




where: 



V(u,v) is a filter function; and 
F(u,v) is said power spectrum. 



A method as recited in Claim 7, comprising: 



determining an angular dependence R6 of said variation R: 



i? 9 (©: * © < © 2 > =7S 2 |F m .J 2 for 0 X s taxr l (-S) < 9 2 
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determining an angular dependence MB of said first moment M: 



M B (Q l Z 0 < d 2 ) --2-2; 



2 2 |F m J 2 



for 0, s tan" 1 (-2) < © 2 . 
1 m 



10 9. A method as recited in Claim 8, comprising: 

dividing said power spectrum into a number of sectors; and 

determining said angular dependence of said variation and said angular dependence of M for each of said 
sectors. 

15 

10. A method as recited in Claim 9, comprising: 

determining a minimum of said angular dependence of M. 

20 11. a method as recited in Claim 4, wherein determining said at least one texture measure comprises: 

determining a root mean square variation of said power spectrum; 
determining an angular dependence of said variation; 
determining a maximum of said angular dependence; 
25 determining a minimum of said angular dependence; 

determining a standard deviation of said variation; and 
determining a relative standard deviation of said variation. 

12. A method as recited in Claim 4, wherein determining said at least one texture measure comprises: 

30 

determining a first moment of said power spectrum; 
determining an angular dependence of said first moment; 
determining a maximum of said angular dependence; 
determining a minimum of said angular dependence; 
35 determining a standard deviation of said first moment; and 

determining a relative standard deviation of said first moment. 

13. A method as recited in Claim 4, comprising: 

40 determining a second power spectrum of a uniform tissue region in said image; and 

normalizing said power spectrum using said second power spectrum. 

14. A method as recited in Claim 13, comprising: 

45 performing edge gradient analysis on said region of interest; and 

determining a maximum cumulative edge gradient. 

15. A method as recited in Claim 4, comprising: 

50 inputting said at least one features into a discriminator; and 

determining a likelihood of risk of fracture using said discriminator. 

16. A method as recited in Claim 11 , comprising: 

55 inputting selected of said variation, said angular dependence of said variation, said maximum of said angular 

dependence, said minimum of said angular dependence, said standard deviation of said variation and said 
relative standard deviation of said variation into a discriminator; and 
determining a likelihood of risk of fracture using said discriminator. 
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17. A method as recited in Claim 12, comprising: 

inputting selected of said first moment, said angular dependence of said first moment, said maximum of said 
angular dependence, said minimum of said angular dependence, said standard deviation of said first moment 
and said relative standard deviation of said first moment into a discriminator; and 
determining a likelihood of risk of fracture using said discriminator. 

18. A method as recited in Claim 1, comprising: 

obtaining a spine image; 

selecting said region of interest in a vertebral body; and 
determining a trabecular pattern in said vertebral body. 

19. A method as recited in Claim 18, comprising: 

background-trend correcting said region of interest. 

20. A method as recited in Claim 18, wherein determining said trabecular pattern comprises: 

determining a power spectrum of said region of interest; and 
determining said at least one texture measure using said power spectrum. 

21. A method as recited in Claim 20, comprising: 

dividing said power spectrum into a number of sectors; and 
determining said at least one texture measure for each of said sectors. 

22. A method as recited in Claim 21 , wherein determining said at least one texture measure comprises: 

determining a root mean square variation and a first moment of said power spectrum. 

23. A method as recited in Claim 1 8, wherein determining said at least one texture measure comprises: 

determining a root mean square variation and a first moment of said power spectrum. 

24. A method as recited in Claim 18, comprising: 

performing edge gradient analysis on said region of interest; and 
determining a maximum cumulative edge gradient. 

25. A method as recited in Claim 18, comprising: 

inputting said at least one features into a discriminator; and 
determining a likelihood of risk of fracture of using said discriminator. 

26. A method as recited in Claim 1 , comprising: 

obtaining a second image of said bone; 

selecting a second region of interest in said second image; 

determining bone mass based on said region of interest and said second region of interest; and 
analyzing said bone based on said bone mass and said texture. 

27. A method as recited in Claim 26, wherein determining said bone mass comprises: 

integrating said region of interest to obtain a first integration value; 

integrating said second region of interest to obtain a second integration value; and 

determining a weighted sum of said first and second integration values. 
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28. A method as recited in Claim 26, comprising: 

obtaining said image and said second image such that each has a plurality of pixels with a gray level value; and 
determining a linear fit from said gray level values of said image and said second image. 

5 

29. A method as recited in Claim 28, comprising: 

determining a first sum of gray levels of first pixels in said image over said region of interest; 
determining a second sum of gray levels of second pixels in said second image over said second region of 
10 interest; and 

determining said bone mass based on said first and second sums and said linear fit. 

30. A method as recited in Claim 26, comprising: 

is obtaining said image at a first energy level and having a plurality of first gray-level pixels; 

obtaining said second image at a second energy level higher than said first energy level and having a plurality 
of second gray-level pixels; 

integrating said image using said region of interest to obtain a first integration value; 
integrating said second image using said second region of interest to obtain a second integration value; 
20 determining a linear fit from said first and second gray-level pixels; 

determining first and second weighted sums of said first and second integration values; and 
determining said bone mass using said first and second weighted sums and said linear fit. 

31. A method as recited in Claim 26, comprising: 

25 

inputting said at least one texture measure and said bone mass into a discriminator; and 
determining a likelihood of risk of fracture using said discriminator. 

32. A method as recited in Claim 31 , wherein: 

30 

determining said at least one texture measure comprises: 

determining a power spectrum of said region of interest; and 

determining a root mean square variation and a first moment of said power spectrum; and 
35 determining said bone mass comprises: 

integrating said region of interest to obtain a first integration value; 

integrating said second region of interest to obtain a second integration value; and 

determining a weighted sum of said first and second integration values. 

40 

33. A method as recited in Claim 1 , comprising: 

selecting said region of interest having a plurality of gray level pixels; and 

determining a surface area of said region of interest at each of a plurality of levels of resolution. 

45 

34. A method as recited in Claim 33, wherein determining said surface area comprises: 

determining a slope using said surface areas and said levels of resolution; 
determining a fractal dimension using said slope. 

50 

35. A method as recited in Claim 34, comprising: 

determining said fractal dimension D = 2 - H, where H is said slope. 

55 36. A method as recited in Claim 33, comprising: 

determining a first surface area of said region of interest based upon said pixels; 
selectively combining said pixels to obtain combined pixels; 
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determining a second surface area based upon said combined pixels; and 
determining a fractal dimension based upon said first and second surface areas. 

37. A method as recited in Claim 36, comprising: 

5 

determining a slope using said first and second surface areas and said levels of resolution; and 
determining said fractal dimension D = 2 - H, where H is said slope. 

38. A system for analyzing bone, comprising: 
a memory (2001); 

an image acquisition device (1000; 2000) connected to said memory; 
a region of interest selection circuit (1002; 2002) connected to said acquisition device; 
a Fourier transform circuit (1005; 2020) connected to said region of interest selection circuit; 
a power spectrum determination circuit (1006, 2030) connected to said Fourier transform circuit; 
a trabecular texture measure circuit (1 007; 2040) connected to said power spectrum determination circuit; and 
a discriminator connected to said trabecular texture measure circuit wherein said texture measure circuit com- 
prises means for determining: 

20 a root mean sq uare variation of said power spectrum H and an angular dependence of said variation ; and/ 

or 

a first moment of said power spectrum H and an angular dependence of said first moment. 

39. A system as recited in Claim 38, comprising: 

25 

a superimposing circuit (1010, 2070) connected to said discriminator; and 
a display (1020; 2080). 

40. A system as recited in Claim 38, comprising: 

30 

a background-trend correction circuit (1004; 2010) connected to said region of interest selection circuit. 

41. A system as recited in Claim 38, wherein said texture measure circuit comprises means for determining a root 
mean square variation and a first moment of said power spectrum. 

35 

42. A system as recited in Claim 41 , wherein said discriminator comprises: 

means for inputting at least one texture measure determined by said texture measure circuit; and 
means for determining a likelihood of risk of fracture using said at least one texture measure. 

40 

43. A system as recited in Claim 38, wherein said texture measure circuit comprises means for determining texture 
measures selected from at least one of: 

a root mean square variation of said power spectrum; 
45 a first angular dependence of said variation; 

a first maximum of said angular dependence; 

a first minimum of said angular dependence; 

a first standard deviation of said variation; 

a first relative standard deviation of said variation; 
50 a first moment of said power spectrum; 

a second angular dependence of said first moment; 

a second maximum of said angular dependence; 

a second minimum of said angular dependence; 

a second standard deviation of said first moment; and 
55 a second relative standard deviation of said first moment. 

44. A system as recited in Claim 43, wherein said discriminator comprises: 
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means for merging said texture measures; and 

means for determining a likelihood of risk of fracture using said texture measures. 
45. A system as recited in Claim 38, further comprising: 

5 

a bone mass entry circuit (2004) connected to said region of interest selection circuit; 

an integration circuit (2006) connected to said bone mass entry circuit; and 

a weighted sum circuit (2008) connected to said integration circuit and said discriminator. 

10 
15 
20 
25 
30 
35 
40 
45 
50 
55 
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Angular dependence of RMS variation: 



Angular dependence of First Moment 
of the Power Spectrum: 





M e (0, £0<O 2 ) = 



I EM 2 



for 0, £ tan 



m n 



The angular dependence of the two measures (IRMS and IFMP) 
is examined by dividing the power spectrum into several 
sectors and performing the summations within each sector. 
From studies f we have found that those with fracture 
elsewhere in the spine exhibit a higher minimum value of 
the angular dependence of IFMP. By taking the minimum valu- 
we force the directionality measure (i.e., perpendicular t« 
the trabeculae) for normal patients since the bone 
trabeculae without osteoporosis is assumed to not show a 
"washed-out" appearance and thus the directionality is 
strong as schematically shown in Figures 6A and B. Since 
the trabeculae are not "washed-out" for normal patients, 
their spatial distribution would contain lower frequency 
structures in a direction perpendicular to the trabeculae. 
Osteoporotic patients would tend to exhibit a more 
isotropic distribution due to the washed-out appearance of 
the trabeculae. Edge gradient analysis on the ROI data can 
also be performed to extract the direction perpendicular t< 
the major trabeculae. The angle exhibiting the largest 
cumulative edge gradient is expected to indicate the 
direction perpendicular to the major trabeculae within the 
ROI. In addition, due to the possibility that quantum 
nottle and x-ray scatter may "hide" the underlying texture 
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